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Abstract. 

We investigate the energy relaxation process produced by thermal baths at 
zero temperature acting on the boundary atoms of chains of classical anharmonic 
oscillators. Time-dependent perturbation theory allows us to obtain an explicit 
solution of the harmonic problem: even in such a simple system nontrivial features 
emerge from the interplay of the different decay rates of Fourier modes. In 
particular, a crossover from an exponential to an inverse-square-root law occurs 
on a time scale proportional to the system size TV. A further crossover back to an 
exponential law is observed only at much longer times (of the order A 3 ). In the 
nonlinear chain, the relaxation process is initially equivalent to the harmonic case 
over a wide time span, as illustrated by simulations of the j3 Fermi— Pasta— Ulam 
model. The distinctive feature is that the second crossover is not observed due 
to the spontaneous appearance of breathers, i.e. space-localized time-periodic 
solutions, that keep a finite residual energy in the lattice. We discuss the 
mechanism yielding such solutions and also explain why it crucially depends on 
the boundary conditions. 
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1. Introduction 

The concept of discrete breathers (DB) has been brought to the foreground by recent 
studies on the dynamical properties of anharmonic lattices |Q . DB are space-localized 
time-periodic solutions, with frequencies lying out of the linear spectrum and they 
have been proved to exist in a wide class of models Q. Although a fairly large 
amount of rigorous and numerical work has been devoted to the identification of the 
conditions for their existence and stability, much less is known about their possible 
relevance for the thermodynamic properties of lattices, both at and out of equilibrium. 
The latter case includes transient dynamics, e.g. relaxation to equipartition [||, and 
stationary heat transport Q| . 

An important example of the role played by DB in determining relaxation 
properties has been provided by Tsironis and Aubry ||. They have shown that a 
characteristic slow relaxation behaviour of the energy emerges when a nonlinear chain 
prepared in a typical thermalized state, corresponding to a given finite temperature, 
is put in contact with a cold bath. This physical setup can be simulated by a damping 
term acting on a number of boundary particles. The system eventually reaches a state 
where a residual finite amount of the initial energy is kept under the form of DB. In 
a subsequent numerical work || it has been further argued that the energy relaxation 
should obey a stretched-exponential law in time. 

Although such a behaviour is considered on its own a signature of complex 
dynamics, its origin remains not completely understood, even in the more widely 
studied problem of glassy relaxation |Q. Loosely speaking, it is usually attributed 
to a multiple-local-minima structure in the energy landscape of the system. This 
would originate metastable states capable of trapping the configuration of the system 
during the relaxation, thus giving rise to a very slow decay. In some cases it is also 
assumed that this behaviour is the result of the competition between two (or among 
many) pure exponential processes |?J. Otherwise, it has been proposed to stem from 
a complex activation mechanism, whereby different degrees of freedom "defreeze" at 
later and later stages. This constitutes a hierarchical scheme, with faster degrees of 
freedom successively constraining the slower ones [|| . 

Beyond the possible relevance of these different scenarios, it is important to have 
a definite quantitative description of the relaxation phenomenon in connection with 
energy localization. With this paper we aim to settle these issues on a firm basis, 
showing also that the picture is more subtle than one might naively expect. As a first 
step, we consider the simple but very instructive case of the harmonic chain, with 
damping acting on its edge particles. This exercise shows indeed how a slower-than- 
cxponential relaxation arises as a global effect from the existence of different time 
scales of the Fourier modes. As a matter of fact, only a linear system with damping 
on each particle would display simple exponential relaxation. More specifically, a 
perturbative calculation reveals that the energy decay law exhibits a first crossover 
from an exponential to an inverse-square-root law. The corresponding crossover time 
To is inversely proportional to the product of the damping constant times the fraction 
of damped particles. 

These results turn out to be extremely useful to be compared with the nonlinear 
case. Remarkably, the only difference concerns the long-time behaviour of the decay 
process. Indeed, in the linear system we observe a second crossover back to an 
exponential law, while in the nonlinear one the energy converges to a finite asymptotic 
value, due to spontaneous localization of the energy in the form of DB. However, we 
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find that finite size effects both in space and time may prevent energy localization, 
making the interpretation of the energy decay process more subtle and interesting. 

The paper is organized as follows. In section || we describe the general features of 
the relaxation dynamics and illustrate the typical layout of numerical simulations. The 
results obtained for the linear system are discussed in section ||; section ^ is devoted 
to the study of the relaxation process in the ft Fermi-Pasta-Ulam (henceforth FPU) 
model (e.g. a quartic nonlinearity in the interparticle potential) and is known to allow 
for the existence of breather states. Finally, we summarize and draw our conclusions 
in section [| 

2. Generalities of the relaxation dynamics 

We consider a homogeneous chain of N atoms of mass m and denote with u p the 
displacement of the p-th. particle from its equilibrium position at time t (for the sake 
of simplicity in what follows we omit the explicit dependence on time). The atoms 
are labelled by the integer space index p — 0, 1, . . . N — 1. In the nearest-neighbour 
approximation of the interactions the equations of motion read 

N-l 

mup = V'(up+i - Up) - V'(u p - Wp-i) - m 2J T pp ,u p >, (1) 

p'=0 

where V(x) is the interparticle potential, for which we assume that V'(0) = 0, denoting 
with V'(x) the derivative of V(x) with respect to x. The last term in eq. (Q) represents 
the interaction of the atoms with a "zero temperature" heat bath in the form of a linear 
damping, characterized by the coupling matrix T. In what follows we consider the case 
in which dissipation acts only on the atoms located at the chain edges, so that T has 
the form 

r pp ' = 7<V,p [^p\o + <V,iv-i] j (2) 

where 7 is the damping constant and o~ p > yP is the Kronccker delta. As we are dealing 
with systems in a finite volume we impose either free-end (u_i = uo, un-i = un) or 
fixed-end (tt_i = un =0) boundary conditions (BC). 

The numerical investigations reported in this paper have been performed 
according to the following procedure. We integrate equations (|l|) starting from a 
microcanonical equilibrium condition corresponding to a given energy density E(0) [N. 
This condition is obtained by letting the hamiltonian system (7 = 0) perform a 
short microcanonical transient, whose duration is typically O(10 3 ) time units. The 
initial condition for such transient is assigned by drawing momenta from a Gaussian 
distribution corresponding to the desired energy, while all displacement variables 
u p are set to zero. During the hamiltonian transient the equations of motion are 
integrated by means of a 5-th order symplectic Runge-Kutta-Nystrom algorithm || , 
while a standard 4-th order Runge-Kutta algorithm is used for the dissipative 
dynamics. 

It is useful to introduce the so-called symmetrized site energies 
e p = -mup + - [V(u p+1 - u p ) + V(u p - tip-i)] . 

The istantaneous total energy of the system is then given by 
N-l 

E=J2 e v ( 3 ) 

p=0 
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Due to the presence of dissipation, the initial energy E(0) decreases in time. Since we 
are interested in determining the decay law of the normalized quantity E(t)/E(0), it 
is convenient to introduce the following indicator 

V(t)=log[-\og(E(t)/E(0))]. (4) 

By plotting T> versus time in a log-lin scale, a streched-exponential law of the form 
E(t)/E(0) = exp[— {t/T) a ] becomes a straight line with slope (3 that intercepts the 
y-axis at — /31og(r). As we shall see later on, this representation is very useful for 
identifying pure exponential regimes (a — 1) and the crossover to slower decay laws. 

Since we are also interested in the effects of DB formation on the relaxation 
process, we introduce a localization parameter C, which provides a rough estimate of 
the degree of energy localization in the system. We define it as 

-,2/ 



C(t)=N—Z -j. (5) 



E 

L P 



e p (t) 



Accordingly, the fewer sites the energy is localized onto, the closer C is to N. On the 
other hand, the more evenly the energy is spread on all the particles, the closer C is 
to a constant of order 1. For instance, in the latter case C for the FPU potential lies 
within the interval [7/4, 19/9], the two extremes being the energy-independent values 
of the harmonic and pure quartic potentials, respectively || . 

In order to smooth out fluctuations, both T>{t) and C(t) have been averaged over 
different realizations (typically 20) of the equilibrium initial conditions. 



3. The harmonic chain 



In this section we discuss the results obtained for the ID harmonic lattice, 
V(x) = ^k 2 x 2 , 

with both free-end and fixed-end BC. In order to solve the problem, we can use 
ordinary time-dependent perturbation theory, provided the damping constant 7 is 
small enough. We note that, since the number of the damped particles is fixed, the 
perturbative approximation is expected to improve by increasing the system size N. 
In the following we shall adopt adimensional units such that k-z = m = a = 1, where 
a is the lattice spacing. As a consequence, time is measured in units of 1/ujq, where 
luq = ^jAk2/m is the maximum frequency of the linear spectrum. 

Let — u)\ and n a (a — 0,1,... N — 1) denote the eigenvalues and normalized 
eigenvectors, respectively, of the unperturbed hamiltonian problem, where 

2 A ■ 2 / la \ 
u a = 4 sm ( y j . 

and the wave-number q a is defined in the Appendix (see eqs. A. 6 and A. 7) . 

In the spirit of time-dependent perturbation theory we look for solutions of the 

form 

JV-l 

u p (t) = c a {t)e-^%. (6) 

a=0 
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where r/™ is the p-th component of the eigenvector n a . By substituting expression (g) 
in the equations of motion, we can rewrite them in the form of a system of differential 
equations for the c a 's. We can then find approximate solutions by expressing the 
latter as power series in the adimensional perturbation parameter j/ujq. The details 
of the calculation are reported in the Appendix. To first order in ^j/luq we obtain 

(7) 



where 



and 



1 1 



(<) = c a (0)exp 

sm 2 (q a ) (9) 

for the free-end and the fixed-end system, respectively, where To — N/2'y. We also 
note that, as a result of the calculation, the damping term introduces a correction of 
the normal frequencies. However, it is straightforward to verify that such a correction 
is of second order in the perturbation parameter and, accordingly, it can be neglected. 

As we see from equations (||) and (^|), the spectra of the decay rates arc 
substantially different for the free-end and fixed-end BC. In the former case, we 
see that the least damped modes are the short-wavelength ones (a ~ N) , the largest 
lifetime being t n _ 1 w 2N 3 /ir 2 j, while the most damped modes are the ones in the 
vicinity of a = 0, with t being the shortest decay time. On the contrary, for 
the fixed-end chain the most damped modes are the ones around a ss N/2, while 
both the short-wavelength and long-wavelength modes are very little damped, being 
T N _! Pa A^(A^+ l) 2 /2vr 2 7 . 

Using the definition of the system energy (||) and the solution (fjj), we can 
cxplicitely evaluate the normalized system energy as 

E2 / n \ 2 -2t/T Q 



E(0) £e£(0) Wl 



(10) 



As we are interested in the typical behaviour of the system when it evolves from 
equilibrium initial conditions, we replace ^(O)^ 2 with its average value 2E(0)/N. 
Thus, recalling equations (Q) and (^|) and approximating for large N the sum over a 
with an integral (r a — > r(g)), we have 

|^ = Ife-^ = e-W^) (11) 



25(0) 

where 2 is the modified zero-order Bessel function. 

Surprisingly enough, although the mode relaxation rates are different in the free- 
end and fixed-end systems, the global result for the energy decay turns out to be the 
same. The asymptotic behaviour of the function ( |TT| ) is given by 

-r t l Ta for t <C r , 



E(t) 



1_ fn, * r . (I 2 ) 
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Hence, r sets the time scale of a crossover from the initial exponential decay, led by 
the longest lifetime, to a power-law decay. 

This peculiar behaviour is illustrated in figure |l]. In order to check the validity of 
the above approximations we compared the analytical result ( |l2] ) with the outcome of 
numerical simulations. We plot in a log-lin scale V vs time for two different lattices 
of size N = 32 with fixed-end and free-end BC. We obtain an excellent agreement at 
least up to t t>_i, where a further crossover to an exponental decay occurs. This 
is an expected finite size effect. Indeed, for times larger than t n _ x only the modes 
around a = N — 1 are populated (in the fixed-end chain also those close to a = 0). 
Therefore, the long-term behavior of the function E(t)/E(0) will be exponential, with 
time constant r„_ 1 /2. However, we note that the second crossover is BC-dependent, 
since the values of t n _ 1 , although both proportional to N 3 , are different in the two 
cases. 

To conclude this section, let us remark that the results reported here apply in their 
present form also to a more general class of models with harmonic on-site potential. 

4. The FPU chain 

In this section we study the FPU model, 

V(x) = \x^ + \x\ (13) 

where the coupling constant of the quartic term has been set to one and we consider 
the energy as the only independent parameter of the hamiltonian system. 

In figure || we plot V vs time for two chains of size N = 256 with free-end and 
fixed-end BC. In both cases no evidence of stretched exponential decay is found. The 
curves follow to a very good extent the characteristic trend of the linear model (see 
equation ( pT| ) ). This remarkable observation can be interpreted by the following 
arguments. First of all, an approximate description in terms of an effective harmonic 
model with energy-dependent renormalized frequencies proved to succesfully account 
for a number of equilibrium properties of the FPU chain |ll| . If we assume that 
the energy extraction rate is slow enough for the system to evolve in a quasi-static 
fashion, we may reasonably expect that a similar scenario as the one described for the 
harmonic case emerges. Second, the harmonic approximation becomes increasingly 
accurate as time elapses, simply because more and more energy is extracted from the 
system by the reservoir. 

Actually, the energy-dependent corrections to the linear spectrum only induce 
in the energy decay curves small deviations from the linear case. Indeed, a closer 
inspection of figure || reveals that changing the initial energy density E(0)/N only 
affects the first crossover region, for both free-end and fixed-end BC. 

A dramatic deviation from the linear-like behaviour emerges however at a later 
stage, when the localization of energy sets in. We observe the spontaneous birth of 
breathers, analogous to what reported in ref This process pins the energy at 
some sites in the chain, thus forcibly causing a further slowing down of the energy 
decay. Nevertheless, we note that this can hardly be detected from the energy decay 
curves, which are mainly determined by the relaxation of the Fourier modes. On 
the other hand, localization is clearly revealed by plotting the parameter C defined 
by formula (|J) (see figure ||). Remarkably, in the free-end chain we have observed 
spontaneous excitation of breathers for all the performed numerical experiments. On 
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the contrary, localization turns out to be strongly inhibited for the fixed-end chain. 
As a matter of fact, we have observed it in a very small fraction (say 5%) of the 
numerical experiments. Moreover, the degree of localization (i.e. the asymptotic value 
of C) depends on the initial energy. 

The pathway to localization is similar to the one described in ref || , whereby one 
breather is formed from a collection of short-lived localized vibrations that eventually 
merge together. The situation is shown for a typical simulation in figure |] (b) 
in the form of a space-time contour plot of the symmetrized site energies. The 
localized modes that emerge out of the relaxation process are very similar to the ones 
that originate from the modulational instability of the zone-boundary modes in the 
hamiltonian system [|[ p^j . This is seen by looking at the displacement and velocity 
pattern of the particles and their spatial spectra. Moreover, the localized objects are 
seen to either stay fixed or move with a definite velocity |fT2f . As a consequence they 
may collide with the boundaries, thereby loosing some energy. However, these losses 
turn out to be tiny, not only because the interaction with the reservoirs are quasi- 
elastic (at least for 7 •< 1) but also because the breather velocity is typically smaller 
than the sound speed. 

To further confirm that the damped system supports breathers similar to the 
ones known from the hamiltonian case, we have performed the following numerical 
experiment. A pattern corresponding to a zero- velocity DB solution in a lattice of size 
N = 256 has been computed numerically within the framework of the rotating wave 
approximation (RWA) This solution has then been used as the initial condition 
for the evolution equations (Q). The results are summarized in figure g, where we 
plot the energy decay curve alongside with a snapshot of the lattice displacements 
at the end of the simulation and the localization parameter. The dissipation acts 
effectively in eliminating the vibrational junk radiating from the approximate solution 
and stabilizing it. The extremely slow decay of the breather energy (see lower panel 
in figu re |q ) is presumably due to anelastic scattering with small amplitude plane 
waves ]14|. 

We can understand what difference the BC make for localization by recalling what 
we have learnt from the harmonic chain together with the known effect of modulational 
instability of the zone-boundary modes jl2|, |l5|. In fact, we see that the main difference 
between the two cases is in the relaxation rates of the Fourier modes (see equation (||) 
and (||)). This implies that in the free-end system the long- wavelength modes quickly 
disappear from the lattice, leaving there just the long-lived zone-boundary ones. This 
explains why modulational instability is so effective for the formation of breathers. 
On the contrary, the "longevity" of the long-wavelength modes in the fixed-end chain 
spoils such instability process, unless a very big energy fluctuation overcomes the other 
effect providing an alternative localization pathway. 

As already pointed out, only the long-term portion of the energy relaxation curves 
is affected by localization. Indeed, in the absence of breathers, the energy slowly drifts 
back to an exponential law. This is clearly seen in figure^, where we plot T> versus time 
for a fixed-end FPU chain. In order to better illustrate this point, we show in figure [| 
(a) two fits to the first power-law portion of the energy decay curves for two values 
of the initial energy. In the lower energy simulation localization has not yet occurred. 
On the contrary, at the higher energy, the same portion of the curve corresponds to 
arousal of localization (see also figure [| (b)). This is reflected in the trend of £ for 
the two cases. We see that the trends of T>(t) are clearly indistinguishable in the two 
cases, making emergence of localization difficult to spot from the energy decay curves. 
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5. Conclusions 

We have analysed the energy relaxation of chains of atoms with linear and FPU 
interaction potential by applying a damping term at the chain boundaries. Both 
free-end and fixed-end conditions have been considered. We have shown that also 
the simple linear system is characterized by nontrivial relaxation features, that are 
determined by the interplay of different relaxation time scales. These are characteristic 
for each Fourier mode, which relaxes exponentially in time with a decay rate depending 
on its wavenumber. In particular, we have derived a two-crossover picture of the decay 
process. More precisely, for t > r = N/2j the initial exponential behaviour slows 
down to an inverse square-root law. At a later stage (t > t n _ 1 oc N 3 /j) the system 
comes back to an exponential decay with a relaxation rate twice the one of the slowest 
Fourier mode. Furthermore, the first crossover is independent of the BC, while the 
second one depends on them. 

We have analysed the nonlinear system under the same conditions and found that 
the behaviour is almost concident with the linear case, provided the energy extraction 
rate is not too fast. Noticeably, our data are not compatible with a stretched- 
exponential relaxation. As a side but important result, we clarify why the fixed- 
end FPU chain does not display spontaneous localization as a product of the energy 
extraction process. This is because the long-wavelength modes of the fixed-end chain 
at equilibrium are characterized by high amplitudes and vanishing relaxation rates, 
that are proportional to -f/N 3 . As a result, the modulational instability of the short- 
wavelength modes is not effective in producing localized vibrations and the energy 
decay behaviour of the fixed-end FPU chain is practically equivalent to the linear 
case. 

The dependence of the characteristic time scales on the system size stems from 
the choice of the damping matrix T. More generally, if one fixes the fraction / of the 
damped particles, the first crossover time tq becomes independent of N (and inversely 
proportional to /), while the second one t n _ 1 is still expected to scale as N 3 . 

In summary, we have brought forward the intrinsic complexity of spontaneous 
breather excitation in spatially extended nonlinear systems. In particular, we have 
illustrated the importance of the system size and of the boundary conditions for such 
studies. 

We want also to observe that the above conclusions hold for ID nonlinear chains, 
where relaxation produces localization of the energy for arbitrarily small initial energy. 
The situation may be significantly different in 2D due to the presence of an energy 
threshold for breather solutions |l6|. Preliminary numerical studies of the 2D FPU 
model seem to confirm such expectations. This point certainly deserves deeper 
investigations, that might provide further insight on the relevance of breathers as 
spontaneous excitations emerging from a relaxation dynamics in nonlinear lattices. 
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Appendix A. 

In this Appendix we report the details of the calculation of the damping rates for the 
harmonic lattice. In matrix notation the equations of motion read 

u = Ku-Tii (A.l) 

where u is the column vector (u , U\, . . . Ujv-i), while K is the matrix 

Kpp' = — [Sp\ P +i + <$p',p-i — 2<5 P ' iP ] (A. 2) 

The top and bottom diagonal elements of matrix K are defined as Kqq = 
if(jv-i),(jv-i) = —k2/m for free-end BC, while definition ( A. 2 ) applies to fixed-end 
BC. The eigenvectors of the hamiltonian problem satisfy the orthonormality conditions 

N-l 

p=0 
N-l 

p,p'— o 



(A.3) 



By substituting expression (|6|) in equation (Al), we get 

N-l N-l N-l 

[c a (t) - u? a c a (t) ~ 2tu; a c Q (t)}e-^ t V ; ~ CaWe"^' K pp< = 

Q=0 Ot=0 p' = 

JV-1 N-l 

- Y, - iwaCcWle-*"-* 51 Tpp'^p" 

By taking the scalar product of the latter equation with n a and recalling the 
orthonormality conditions contained in ( |A.3| ), we are led to the system of equations 

N-l 

c a (t) - 2iu a b a {t) = ~Y [Cc'(*) - »Wa'C a /(t)] r^'e^^'-^)*, (A.4) 

a'=0 

where 

r QQ ' = (r, Q ,r^')= E VpTpp'^ (A.5) 

p 7 p'— o 



We can obtain approximate solutions of system (A.4) by means of a perturbative series 
of the form 

where ct if), Ca (t), . . . are first order, second order amplitudes, and so on, in the 
perturbation parameter 7. We can then use the usual iterative procedure. If we assume 
that at t = only the mode a = a is populated, we can approximate c a (t) = 5 a . & 
(independent of t) in the right-hand side of equation (A.4) and then integrate to get 
Ca\t). The procedure can be repeated for obtaining (£), etc. For the first order 
approximation we obtain the following equation 

c Q (f) - [2icu a - T aa ]c a {t) - iuj a T aa c a (t) = 0, 

which is readily integrated, yielding expression (^). 
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To calculate the matrix elements T aa , we need the explicit expression of the 
eigenvectors of the unperturbed problem, which are 



n. 



p 



p V N 
for the free-end chain, and 

Vp = v/^ sm [ga(p+ 1)] , 



3q 



N 



q a = 



(a + l)ir 
N+ 1 ' 



a = 0, 1...N-1 



a = 0A...N -I 



(A.6) 



(A.7) 



for the fixed-end one. By inserting these expressions in equation (A. 5), one can easily 
obtain equation (N) and equation (ph. 
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(a) (b) 



Figure 1. Symbols represent T> versus time for a harmonic chain with 
N = 32, and with free-end (a) and fixed-end (b) boundary conditions, ■y = 
0.05,0.1, respectively. Solid lines are plots of formula ^L~\) and of the function 

Cxp(-2t/"TiV-l). 
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Figure 2. Symbols represent T> versus time for an FPU chain with N = 256, 
7 = 0.05 and with free-end (a) and fixed-end (b) boundary conditions. Solid line 
is the plot of formula nlA). 




Figure 3. The localization parameter C vs time for an FPU chain of size 
N = 256 with free-end (a) and fixed-end (b) boundary conditions. 
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Figure 4. FPU chain with N = 256. (a) Lower panel: log-log plot of the 
first power-law portion of the energy decay curves. Best fits (solid lines) give 
exponents —0.51 and —0.53 for E(t)/E(0) = 0.05 and 1.0, respectively. The 
upper panels show the trend of the localization parameter C for both cases in the 
same time domain, (b) Space-time contour plot of the symmetrized site energies 
for E(t)/E(Q) = 1.0. The figure shows a close-up in the region of the chain where 
the breather develops. 
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Figure 5. Plots of an RWA solution relaxing in an N = 256 chain. Lower 
panel: normalized energy vs time. Upper left panel: snapshot of the atomic 
displacements at t = 5 ■ 10 6 . Upper right panel: the localization parameter C vs 
time. 
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Figure 6. Symbol represent T> versus time for an FPU chain with N = 16 and 
fixed-end boundary conditions. Solid and dashed lines are plots of formula (Wjlj 
and of the law exp(— 2t/rjv— i ), with rjy-i = N(N + l) 2 /27r 2 7, respectively. 



